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Kalman filters have been used for some time to process clock measurement 
data. Due to instabilities in the standard Kalman filter algorithms, the results 
have been unreliable and difficult to obtain. Often, in order to obtain reasonable 
results, the data has had to be manually edited, the filter fine tuned, or the model 
adjusted by the analyst. During the past several years, stable forms of the Kalman 
filter have been developed, implemented, and used in many diverse applications. 

These algorithms, while algebraically equivalent to the standard Kalman filter, ex- 
hibit excellent numerical properties. Two of these stable algorithms, the Upper 
triangular-Diagonal (UD) filter and the Square Root Information Filter (SRIF), 
have been implemented to replace the standard Kalman filter used to process data 
from the DSN hydrogen maser clocks. The data are time offsets between the clocks 
in the DSN, the timescale at the National Institute of Standards and Technology 
(NIST), and two geographically intermediate clocks. The measurements are made 
by using the GPS navigation satellites in mutual view between clocks. The filter 
programs allow the user to easily modify the clock models, the GPS satellite de- 
pendent biases, and the random noise levels in order to compare different modeling 
assumptions. 

The results of this study show the usefulness of such software for processing 
clock data. The UD filter is indeed a stable, efficient, and flexible method for 
obtaining optimal estimates of clock offsets , offset rates, and drift rates. A brief 
overview of the UD filter is also given. 

I. Introduction ordinated Universal Time (UTC) as realized at the Na- 

tional Institute of Standards and Technology (NIST), with 
The requirements of the DSN-complex clocks are that a knowledge of 1 microsecond and that the rate of each 
each clock be maintained within 10 microseconds of Co- clock be kept within ±IE - 12 df/f with a knowledge of 
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±3E — 13 df/f, The quantity df/f is known as the frac- 
tional frequency deviation where 


df/f = 


and f m is the measured frequency and /o is the reference 
frequency. Also there is a requirement that a permanent 
record be kept of the synchronization and syntonization 
of the DSN-complex clocks. In order to meet these re- 
quirements, the time and frequency offset to UTC (NIST) 
is monitored using the Global Positioning System (GPS) 
navigation space vehicles. The goal (not yet realized) is 
to adjust the rates of all the clocks at the same time at 
approximately 200-day intervals. The clock rates (frequen- 
cies) are adjusted when they are above UTC (NBS) by a 
certain amount to below, and then allowed to drift through 
zero rate offset from UTC (NIST) until they are again high 
enough to require adjustment. This technique, done with 
some care, allows one to keep the clocks within the re- 
quired time and frequency offset limits and still not have 
any jumps in the time offset. All of the drift rates of 
the hydrogen masers are about the same, therefore, the 
complex clocks can be kept in close synchronization and 
syntonization. 

Monitoring of the DSN-complex clocks is done by the 
operations analyst at JPL. Adjustments of the clocks is 
recommended by the operations analyst with the advice 
of the engineering staff. The adjustment is scheduled af- 
ter a check of the users of the DSN and is made during 
a maintenance period at the complex. The actual adjust- 
ment is made by the engineering staff at the complex. It 
is assumed that the drift of the hydrogen maser rate is 
due to cavity drift, therefore, rate adjustments are made 
by adjusting the varactor, thereby changing the apparent 
cavity size of the hydrogen maser. 

At critical times during some projects (such as the 
Voyager Uranus encounter), it is desirable to have the 
clocks synchronized and syntonized more closely than the 
above specifications. When this is the situation, the oper- 
ations analyst will predict the time and rate of the DSN 
clocks 30 to 90 days prior to the encounter date and set 
the rate to allow the clocks to drift to the near zero time 
and rate offset by the critical time. At the time of the 
Voyager Uranus encounter, JPL engineering staff used a 
Kalman filter to predict the clock offset 60 days into the 
future with excellent results [3]. 


At JPL, a Kalman filter has been used to process 
clock offset measurements but it has never been used oper- 
ationally. The present operational environment for clock 
management in the DSN is a personal computer with a 
considerable amount of manual intervention. The data is 
handled in weekly batches. By 1992 the process is to be 
transferred to the Network Frequency and Timing subsys- 
tem (NFT), which is a part of the DSN Frequency and 
Timing system (DFT). The flow of data and estimate up- 
dates will occur more often (perhaps as often as once per 
hour) so as to allow operations to be able to monitor clock 
performance in near real time. This will require a fairly 
robust processing environment with automatic elimination 
of outliers and other techniques to make the process as 
autonomous as possible. It is expected that operations- 
analyst intervention will be relegated to filter adjustments 
in response to scheduled changes of time or rate in any of 
the clocks being entered into the filter parameters. 

The operations analysts will have the opportunity to 
do postanalysis of the data. In the near future, it is ex- 
pected that some data will continue to be available on a 
weekly basis and that better ephemerides will be available 
several weeks after the fact. At this time, use of the fil- 
ter/smoother by the operations analyst can produce better 
estimates of the clock parameters for archiving. To achieve 
these goals, a robust filter and smoother program is under 
development at JPL. The results reported in this article 
were obtained by running this program with data as de- 
scribed in the next section. No “fine tuning” adjustments 
or manual data editing were performed. The output, how- 
ever, was excellent. The filter used in the JPL program 
is the numerically stable Upper triangular-Diagonal (UD) 
form of the Kalman filter. The UD filter has been under 
development for several years [1] and has earned a repu- 
tation for being accurate and efficient. The more tradi- 
tional Kalman formulation has been known for some time 
to be numerically unstable. That is, due to rounding errors 
within the computer, the results produced by the filter pro- 
gram are sometimes completely wrong. The smoother used 
in this code is a recent improvement [2] of the Rauch-Tung- 
Streibel (RTS) smoother [5]. The revised RTS smoother 
is also robust and accurate. 

The filter and smoother programs make use of subrou- 
tines from the Estimation Subroutine Library (ESL). This 
is a collection of FORTRAN77 subroutines for construct- 
ing filters and smoothers. The routines in the ESL have 
been carefully coded and tested to ensure their accuracy 
and reliability. Overall, the goal is to produce optimal es- 
timates of clock parameters without manual intervention. 
The results of this first attempt are extremely encouraging. 


191 


II. Clock Measurements 


There is an NIST-designed GPS timing receiver lo- 
cated at each complex. A timing pulse from the complex 
master clock is fed to each receiver so that the data output 
from the receiver is the offset of the complex master clock 
to GPS timescale. The receivers are maintained on the mu- 
tual view schedule, which is generated by the Bureau Inter- 
national des Poids et Measures (BIPM) in Paris, France. 
These data are kept in a database along with data from 
UTC (RRL) and UTC (TAO) in Japan, UTC (NIST), 
and the NIST cesium clock at WWVH Hawaii, U. S. The 
data from Japan are obtained from the GE MKIII cata- 
logue, which is administered worldwide by the BIPM and 
inside the U. S. by the United States Naval Observatory 
(USNO). Each of the data lines is tagged to indicate the 
receiver from which it came and then placed in a file and 
stored on a disk in a personal computer. There is one file 
for each day, and it is named by Modified Julian Day. 1 

The first process is to take the first difference of the 
data to obtain the differences between the clocks on the 
ground, thus eliminating the space vehicle clocks. There 
are six clocks used; because of the geometry of the clock 
locations there are 11 possible mutual views available 
(Fig. 1). A program is used to produce a file which con- 
tains as entries the station pair, the satellite number, the 
MJD, the time (second) the data was collected, and the 
difference in time in tenths of nanoseconds between the 
first and second clock. (The above format could easily be 
adapted to include clock difference data from other sources 
such as two-way time coordination.) Data was generated 
by this program from MJD 47000 to 47132. There were 
approximately 4000 measurements in the data set. 


d{(t): the drift rate, that is, the second derivative of 
the offset, d, — d 2 Si/dt 2 . 

&*w(0 : the measurement biases. These depend on the 
two stations, i and it, and on the satellite, /, being ob- 
served. These biases will be explained in more detail 
• later. 

The dynamics model describes how the above quan- 
tities change over time. This is expressed as follows: 


S|(0 = r *(<) 

(l) 

r,(0 = M 0 

(2) 

•to 

tl 

*"13 

(3) 


where • denotes time differentiation, and tu,*($) is white 
noise, with zero mean and known power spectral density 
(psd) W . That is, 


E 


i^(f)u;i(r)J = W6(t — r) 


(4) 


where E [ •] is the expected value operation, and where 6 is 
the Dirac delta function. Note that Eqs. (1) and (2) are 
just restatements of the definitions of r, and d,-. Equa- 
tion (3) says that the rates are “approximately” constant. 

In order to apply a discrete filter to this problem, 
Eqs. (1-3) must be put into the form of a discrete dynamics 
equation 


III. Mathematical Model 

In order to describe the mathematical model used for 
this analysis, the following notation is introduced: 

Si(t): the clock offset at station i. That is, $,• = c,* — T, 
where c, is the clock reading at station i and T is the 
“true” time. For this analysis, “true” time is taken to 
be the time given by UTC (NIST). 

ri(t): the clock offset rate, also called the frequency 
error. This is just the time derivative of the offset, 
i.e., r< = dsi/dt. 

1 The Modified Julian Day (MJD) is a continuous day count 
with an initial epoch of 0000 hours UT on November 17, 1858. 


x(t + At) = $(t, At)x(t) + tu(t) (5) 

where x denotes a vector with components that include 
s,-, **«, di for each station. Note that all of these quanti- 
ties are time-dependent; however, to simplify the notation, 
this dependence will often be suppressed. The matrix 4>, 
called the transition matrix, is determined by writing out 
a discrete form of Eqs. (1-3). A simple way to do this, 
for example, at least for the offsets and rates, is to use a 
discrete form of Taylor’s Theorem: 

Si(t + A0 = Si(t) + n(t)At + 0M(t)At 2 (6) 

r,(t + AO = r,(0 + di(t)At (7) 
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For the drift rates, this same approach produces 

/■A< 

di(t + At) = d{(i) + J Wi(r)dT 


( 8 ) 


where Wi(t) represents the random noise. However, the 
effects of the noise on the drift rates are shown directly 
in Eq. (8), but this same noise also affects the offsets and 
frequency errors because of the presence of d{(t) in Eqs. (6) 
and (7). As written above, however, there will be a delay, 
caused by the time discretization, before the noise at time 
t propagates into the offset and rate. To avoid this delay, 
noise terms are included directly in the equations for Si 
and r,-. This is done by using the continuous form of the 
transition matrix to propagate the noise over the interval 
t to t -f A* . That is, if the noise that enters the continuous 
system at time t is given by the vector 


w(t) = (o,0,uv(0) 

where u; r (f) is white noise with zero mean and psd W r} 
then the integrated effect of the noise at time t + At will 
be the vector 

q(t + A<) = ( q„q r ,qd) T 

rp 

where ( q ,, q r ,q t ) ~ N(0,Q) and 
Q = E $ (At — t) w(r)dr J w T (<r)$ T (At — a)do 

= E [fL 


$(r)w(At — t)w j (At — a)# 1 (a)d<rdT 


= (0.5r 2 , r, l) T W r (0.5r 2 , r, l)dr 


( 9 ) 


since 


*(<) = 


1 t 0.5 t 2 
0 1 t 

0 0 1 


and E[wJ w; r ] = W r . The integration in Eq. (9) gives 


Q=w r 


r^A< 5 §A t 4 I At 3 
|A t 4 |A t 3 \At 2 


L t 3 2 At 2 


At 


( 10 ) 


Thus the discrete dynamics model for these filter states is 


where 


/ Si(t + At)) 


/«.(<) \ 

r,(< + A<) 

= 

»•.'(<) 

\ d{(t + At) J 


\di(t)J 


'1 At 

*A< 2 ‘ 

= 

0 1 

At 


.0 0 

1 


+ u>i 


(ii) 


( 12 ) 


Wi ~ N(0 , Q t ), and Qi is given by Eq. (10). There will be a 
transition, of the form of Eq. (11), for each station. Thus, 
the full transition matrix $ in Eq. (5) will have blocks of 
the form of Eq. (12) on its diagonal. The process noise 
term, w(t) i in Eq. (5) will have as its covariance matrix a 
block diagonal matrix, with 3x3 matrices of the form of 
Eq. (10) on the diagonal. 

The measurements to be used for estimating the off- 
sets, drifts, and drift rates are the differences in clock val- 
ues from the GPS times, as measured at two different sta- 
tions. That is, the measurements are 


z ikl (f j ) 


which are assumed to be of the form 


*.«(<;) = (ci(tj) - GPS,) - - GPS,) 

+ &»*/(*>) + Vik(tj) 

= (f J ) ” c k(tj) bikl(tj) H" v ik(tj) (13) 


where, as defined earlier, c,(fj) is the clock reading at sta- 
tion i at time tj. The fact that s;(i,) = c,(tfj) — NIST(tj) 
allows one also to write this as 
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Zik,(tj) = - NIST(t,)) - (c k (tj) - NIST(tj ) ) 

+ bikl(tj ) + v ik(tj) 

= s k(tj} "I" 

The first expression is used to obtain the measurement 
data as described in Section II. That is, the data actu- 
ally used for the analysis are values for c,(< ; ) — GPSj(tj). 
However, for the filter/smoother, the measurements are as- 
sumed to be of the second form, where the biases, 6,*/(tj), 
satisfy the trivial dynamics equation 

b ikl (t + At) = b ikl (t) (14) 

and the measurement noise, Vi k (tj), is, as usual, a nor- 
mally distributed zero mean random variable with con- 
stant variance r. The biases are appended to the “core 
states” to obtain the actual state vector used in the filter: 

x = (si, ri, di, S2, r2, c?2, ...» Sm, r m ,d m , 6ni, 6112, 6113, . . .) 

The transition matrix will have an identity matrix at- 
tached to it to reflect Eq. (14). Then the measurements 
can be written in standard form as 

Zj = HjXj + Vj (15) 

where Hj is a 1 x n measurement matrix, consisting, in 
this case, of two positive ones and a negative one, in ap- 
propriate positions. For example, a measurement based 
on observations from stations 1 and 3, and satellite 6, will 
have a measurement matrix 


H = (l i 0,0,0, 0,0, -1,0,0,. ..,0,1,0,...) 

where the final one appears in the position corresponding 
to bias 6 x 36 . 


IV. The UD Form of the Kalman Filter 

The Kalman Filter has become the standard tool for 
computing optimal estimates of a state vector, x(t), gov- 
erned by a mathematical model of the form of Eqs. (5) and 
(15). To describe the form of the filter that is of interest 
here, the following notation will be used: 


x k \j: estimate of x(t k ) using data up to time tj 

P k \j : covariance matrix for the error in x k y as an ap- 
proximation to x(tk) 

The Kalman Filter time update, also called the prediction 
step is 


£jbl*-l — ®kxk-l\k-l 

(16) 

Pk\k-i — $kPk-i\k-i$I + Q* 

(17) 


The measurement update is accomplished via the following 
set of equations; 


Bk — H k Pk\k-\Hk + Rk 


(18) 

* 

II 

>p 

1 

to 
*■ 1 

(Kalman gain) 

(19) 

P k \ k = {I-K h H h )P k 


( 20 ) 

Vk = Zk — HkXk\k-\ 

(innovations) 

( 21 ) 

£jb|jfc = £fc|ib-l + Kk^k 


( 22 ) 


The numerical instability of the Kalman equations, as 
noted in Section I, is due to the structure of these equa- 
tions. The most common indication that something has 
gone wrong is when rounding error causes P k \ k to be not 
positive definite. Note that it is not clear from these equa- 
tions that P k is even symmetric. 

The UD form of the Kalman Filter is algebraically 
equivalent to the above equations, but has much better nu- 
merical properties. It is based on the fact that any positive 
definite matrix P can be written as a product of matrices 

P = UDU t (23) 

where U is an upper triangular matrix with ones on the di- 
agonal, and D is a diagonal matrix, with positive diagonal 
elements. The idea behind the UD formulation is to start 
with the UD factors of P 0 1 0 , then compute the factors of 
all the remaining covariances, without ever computing the 
covariances themselves. It can be shown that by avoiding 
the covariances, the accuracy of the results is greatly im- 
proved. On the other hand, if the covariances are needed 
for output they can be easily obtained by simply multiply- 
ing the factors together. 
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To describe the time update in UD form, assume that 
the factors of P k j* are known from the previous step. That 
is, JJ k \ k and D k \ k have already been computed, where 

Pk\k = U kik D klk U^ k (24) 

Then x k+l \ k is computed just as in the Kalman formula- 
tion; but the equation for P k +i\ k becomes 

R fc+l|Jb = $kPk\k$7c + Qk 

= $kU k \kD k \ k U k \ k <bJ + Q k (25) 

Since $ k is upper triangular, so is the product $kUk\k* 
Hence, if Q k were zero, then the factors of P*+i|jb would 
be just 

RkD k R k 

where 

Rk = *uU klk (26) 

In order to account for the process noise, Q k is factored 

into its UD factors 

Qk = UqDqUj 

which can also be written as 

Qk = Y,d iqi qJ (27) 

1 = 1 

where qi is the ith column of U q . Then Eq. (25), using 
Eqs. (26) and (27) becomes 

n 

Rk+I\k = RkDkRk + y ^diqiqj (28) 

»=1 

Now the time update step is completed by finding the UD 
factors of P Jk _|_ 1 |j fc . This is done by a series of “rank one” 
adjustments to previous factors. That is, consider the gen- 
eral problem: given factors i?, jD, compute the factors of 
the matrix 

RDR T + dqq T 

where d is a positive scalar, and q is a vector. 


The matrix dqq T is a matrix of rank one. Hence dqq T 
is called a “rank one” adjustment to RDR T . There is a 
very accurate and efficient algorithm for computing the 
factors of this new matrix, given the factors R , D , and the 
quantities d and q. It is important here that the factors 
are not multiplied together, as this would destroy the nu- 
merical accuracy. This rank-one update algorithm, called 
the Turner- Agee Algorithm [1], is applied repeatedly to 
Eq. (28) to complete the time update step of the filter. 

Next, consider the measurement update, as given in 
the Kalman form by Eqs. (18-22). The Kalman equation 
can be used to process vector measurements. However, the 
UD form is more efficient if the measurements are scalars. 
In this case, in the measurement update equations, B k is 
a scalar, as is v ki and K is a vector. The measurement 
update is done several times, once for each scalar mea- 
surement for each time update. After the previous time 
update, the factors U k \ k -i and D k \ k _\ of P k \k-i are avail- 
able. The (scalar) B k as defined by Eq. (18) is computed 
as follows: 

B k = H k U k \ k ^iD k \ k ^iU k ^ k ^ 1 H k + Rk 
= yl^k\k-iyk + Rk 

where 

y* = 

is a vector. Next, Eqs. (19) and (20) are written as 
P k]k = (I-I< k H k )P k \ k _ x 

= - j^Rk\k-\Hk HkPk\k-x 

- U k \ k -\D k \ k .iU^ k _ l 

- ~^Uk\k~i^k\k-\U k \k~\Rk H k U k \ k ^\D k ^ k ^iU k ^ k _ 1 

= u k \k-x [am*— i - j- k vyT ] ^*-1 ( 20 ) 

where V is the vector defined by 

V = D k \k-xUlk-xHk 
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All that remains is to compute the factors U v , D v of 

Dk\k-i ~ j-VV T (30) 

Then it follows that 

Pk\k = U k \ k - 1 U v D v UZU? lk _ 1 (31) 

^*|k ^k|k uT.. 


Notice that Eq. 30 is a rank-one modification of a diagonal 
matrix. Hence its factors can be computed accurately and 
efficiently, just as was done in the time update. 

While this UD formulation may seem more compli- 
cated than the relatively simple matrix equations of the 
Kalman form (Eqs. 16-22), a computer program to imple- 
ment the UD filter is very straightforward. Two special 
subroutines are needed: 

(1) A subroutine to efficiently compute the product of two 
upper triangular matrices, such as 

(2) A subroutine to adjust the factors of a matrix by 
a rank-one addition: given factors £/, D } a positive 
scalar c, and a vector y, compute the factors of U DU T 
+ cyy T 

Then the time update is done by the following steps: 

(1) Compute R = $kU (sometimes called the determin- 
istic update) 

(2) For each process noise term d^qj, update the factors 
using the rank-one adjustment routine 

The measurement update is also broken down into a series 
of simple steps, one for each scalar measurement: 

(1) Compute the scalar B = y 7 ^ Dy + r where y = U T H T 

(2) Use the rank-one adjustment routine to compute the 
factors of D - £ VV T where V = DU T H T 

(3) Combine these factors with the previous factors, ac- 
cording to Eq. 31 


V. Filter Implementation 

In addition to the clock measurement data, as de- 
scribed in Section II, the filter routine must be given the 
following information: 

(1) Prior estimates of the state vector 

(2) Uncertainties associated with the state estimates 

(3) Process-noise standard deviations 

(4) Measurement-noise standard deviations 

Generally, the filter is somewhat insensitive to prior 
estimates of the state; hence, they are usually initialized at 
zero. The uncertainties must be large enough so that the 
filter will accept the measurement data and use it to adjust 
the state estimates, except that the offsets for the NIST 
station should remain close to zero. The process noise 
can depend on the station; however, for the results of Sec- 
tion VI, all stations were given the same amount of noise. 
Similarly, the measurement noise can be made a function 
of the stations and the space vehicle. For simplicity during 
this first analysis, however, the measurement noise uncer- 
tainty was assumed to be the same for all measurements. 

In order to run the smoother, certain information 
must be saved at each measurement update time. Thus, 
to limit the amount of storage required to run the filter 
and smoother, the program was set up to process mea- 
surements in so-called “mini-batches.” Each mini-batch 
consists of all measurements obtained during a fixed time- 
span, say one day, or a half a day. The midpoint of this 
timespan is chosen as the time of the update (also called 
the “epoch time”). The measurements within the mini- 
batch are translated to the epoch time by using the same 
transition matrix, $(ti,f 2 ), as is used in the dynamics 
equation. More precisely, suppose that the actual mea- 
surement is 

z(ti) = Hx(ti) + tv (32) 

where t\ denotes the time at which the measurement was 
taken. Let 1 2 be the epoch time for this mini-batch. Then, 
if it is assumed that the timespan for this mini-batch is 
sufficiently short, so that the integrated effect of the noise 
is negligible, one can write 

*(<i) = $(<2 ,<i)x(* 2 ) (33) 

so that Eq. (32) becomes 

z(*i) = H$(t 2 ,ti)x(t 2 ) + w (34) 
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That is, the measurement z(t\) becomes a measure- 
ment of the state at the epoch time, but with the measure- 
ment matrix given by H$ instead of just H. It is impor- 
tant to realize that this use of mini-batches can change the 
filter results because it assumes that the effect of the pro- 
cess noise over the timespan of the mini-batch is negligible. 
The code that was developed allows the user to select one 
or more mini-batches per day, or to process each measure- 
ment as it is received. 

The filter/smoother was implemented in FORTRAN- 
77, making use of routines from the Estimation Subroutine 
Library (ESL) to do most of the numerical computations. 
The UD form of the Kalman filter, as described in Sec- 
tion IV was used for the filter; the smoother made use 
of Bierman’s modification of the RTS smoother [2]. The 
general structure of the program is as follows: 

(1) Read a namelist file to initialize the filter: 

Start and stop times 

Maximum number of measurements to process 
Initial error covariances and state estimate 

(2) For each measurement: 

If new minibatch, 

(a) Do time update 

(b) Print estimates and estimate error standard de- 
viations 

(c) Translate measurements to epoch time and do 
the measurement update; write the smoother 
gains to a file for later use 

(3) After all measurements have been processed, print out 
the correlation coefficients 

(4) If the smoother has been requested, run the smoother 
back to the start time 

The results of Section VI were obtained by running 
the program with 6 stations, 11 station pairs, and 6 satel- 
lites. Hence the filter state vector has 84 components. An 
important feature of the filter implementation is the data 
editing capability. Before each measurement is actually 
included into the estimate, a test of the measurements 
consistency with past measurements is performed. This 
is done as follows. First the residual (or innovations) as 
defined by Eq. (21) is computed, as is the variance of this 
residual which is just the quantity Bk defined by Eq. (18). 


If the square of the residual divided by its variance is less 
than a user-specified tolerance, then the measurement is 
accepted and used to update the current estimates. For 
the results described in Section VI, a tolerance of 400 was 
used. This means that a measurement was rejected only 
if its residual was more than 20 standard deviations away 
from its mean value of zero. 

VI. Numerical Results 

The numerical results described in this section were 
obtained by running the filter and smoother on the data 
described in Section II. Approximately 10 data points were 
rejected by the data editor that is built into the filter. Fig- 
ure 2 gives an overview of the results of the smoother. This 
figure shows the offsets from UTC (NIST) of the three 
clocks in the DSN. This is a useful display for the oper- 
ations analyst. It can be used to tell when adjustments 
to the rates are in order. Notice that a rate adjustment 
would have been appropriate around MJD 47080; however 
an adjustment was not made until MJD 47130. The Allan 
variances of the smoother results, shown in Fig. 3, are con- 
sistent with those reported in [4], but are somewhat better. 
It is felt that these variances can be further improved with 
better modeling. 

In order to contrast the output from the smoother 
with the input to the filter, see Fig. 4(a). This shows 
the raw data for the California clock, with the smoother 
results superimposed. A rate of 1.2 E — 13 was removed 
from both plots. The smoother nicely bridges a gap in 
the data at about MJD 47050. During this gap, the stan- 
dard deviations computed by the smoother increased from 
about 4 nanoseconds, on MJD 47046, to 40 nanoseconds 
on MJD 47052, then dropped back to 4.4 nanoseconds on 
MJD 47053. Other data gaps were also handled well. This 
ability of the filter /smoot her to bridge data gaps will be 
important for operational uses. Figure 4(b) shows a super- 
position of the results produced by the JPL smoother and 
the NIST filtered estimates. Again, the 1.2E—13 rate has 
been removed from both results. The agreement is quite 
good, with less than 8 nanoseconds difference. The con- 
sistent bias, with the JPL results slightly lower, may be 
a result of the way in which the space vehicle biases are 
handled. 

The offset rate of the complex clocks is the most im- 
portant parameter from the operational point of view. 
Presently, the rate is determined by taking the mean of all 
the GPS measurements for ten days, then doing a least- 
squares linear fit and using the slope of that fit as an esti- 
mate of the rate. 


197 



The results of these “hand calculations” are shown as 
dots in Fig. 5. Generally, there is good agreement between 
these results and those produced by the filter/smoother. A 
careful examination of the times at which the hand calcu- 
lations differ noticeably from the smoother results revealed 
two reasons for these discrepancies. For example, in the 
vicinity of MJD 47060-47070, there was very little data 
for the clock in Spain. This caused the smoother to re- 
port a large uncertainty in its rates over this interval, and 
explains the difference between the smoother output and 
the hand calculation shown in Fig 5(c). The point ob- 
tained by the hand calculations for the clock in Australia 
on MJD 47070 was found, on closer examination, to have 
been calculated incorrectly. The correctly calculated value 
is very close to that produced by the smoother. 

VII. Conclusions and Future Work 

The UD form of the Kalman Filter appears to solve 
many of the problems encountered in the past with pro- 
cessing clock data on a routine basis. The numerical stabil- 
ity of the filter/smoother provides an accurate and reliable 
method for computing optimal estimates of the offsets, off- 
set rates, and drift rates. The automatic data editing fea- 
ture removes bad data points appropriately. 


The first several days worth of the data that were used 
for this study seemed to be very poor. For this reason, 
data corresponding to the first four days were skipped. 
As a comparison, this same data was also processed by 
using a Square Root Information Filter (SRIF). The SRIF 
is an even more numerically reliable form of the Kalman 
Filter. This increased reliability, however, is attained at 
the cost of an increase in computational time. The SRIF 
was able to handle the bad data at the beginning of the 
data set with no problems. Since computer run-time is of 
little concern for this application, it may be advisable to 
replace the UD filter with a SRIF. 

Possible improvements to the mathematical model 
used in the filter would be of interest. For example, it 
is clear that the assumption about the biases being con- 
stant is incorrect. An improved model would allow the 
biases to be time varying. The rate at which they should 
vary, however, would have to be determined by statistical 
parameter estimation techniques, such as maximum likeli- 
hood estimation. Also, in the current model, all clocks are 
treated as if they were identical. In fact there are three 
different types of clocks in the DSN. Special characteristics 
of the clocks could be taken into account to improve the 
accuracy of the model. 
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Fig. 5. Kalman smoother output and hand calculation of 
rate using a linear least-squares fit over 10 days of the 
time offset: (a) California, (b) Australia, and (c) Spain. 
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